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ABSTRACT 

We numerically investigate the long-time evolution of density perturbations 
after the first appearance of caustics in an expanding cosmological model with 
one-dimensional 'single-wave' initial conditions. Focussing on the time-intervals 
of caustic appearances and the spatial distribution of caustics at subsequent 
times, we find that the time-intervals of caustic appearances approach a constant, 
i.e., their time-subsequent ratio converges to 1; it is also found that the spatial 
distribution of caustics at a given time features some universality rules, e.g., the 
ratio between the position of the nearest caustic from the center and that of the 
second nearest caustic from the center approaches a constant. Furthermore we 
find some rules for the mass distribution for each caustic. Using these universality 
constants we are in the position to predict the spatial distribution of caustics at 
an arbitrary time in order to give an estimate for the power spectral index in the 
fully-developed non-dissipative turbulent ('virialized') regime. 

Subject headings: cosmology: theory - large-scale structures - self-similar - non-dissipative 
gravitational turbulence 



1. INTRODUCTION 

The evolution of large-scale structure in an expanding universe model is one of the 
most important and interesting problems in cosmology. It is generally believed that these 
structures have formed owing to gravitational instability. After structure has formed, the 
combined action of multi-stream forces and gravity is responsible for their stability in time 
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(Buchert & Dommguez 1998), eventually leading to 'virialized' structures on certain spatial 
scales such as galaxies and clusters of galaxies. Such equilibrium states, if they persist in time, 
are not necessarily in thermodynamical equilibrium, but are the consequence of a process that 
was termed by Gurevich & Zybin (see their review 1995 and refs. therein) non-dissipative 
gravitational turbulence, emphasizing the role of multi-stream forces that counteract the 
gravitational action and so eventually lead to equilibrium structures, supported by velocity 
dispersion in collisionless systems. 

Here we consider collisionless particles as candidates for cold dark matter, and focus our 
interest on universal behavior in the gravitational collapse process. The hope is to discover 
some rules that govern the multi-streaming hierarchy in a self-gravitating system. In this 
research note we shall concentrate on the simplest model of a plane-symmetric 'single-wave' 
perturbation on an expanding universe model; the discovery of universality laws needs a 
detailed study at high-spatial resolution. An advantage of plane symmetry is the existence 
of exact solutions that we shall incorporate in our numerical scheme. Other aspects of 
multi-streaming have already been studied in three spatial dimensions (see, e.g., Kofman 
et al. 1994), and there is a plethora of works on powerlaw, fractal, or self-similar scaling 
properties of the density field, which are all indirectly related to the understanding of multi- 
stream systems. Also numerical simulations of self-gravitating systems practically evolve a 
multi-stream flow on smaller spatial scales. 

Recently, Yano & Gouda (1998) have studied the time evolution of density perturba- 
tions in the plane-symmetric system. They showed that the power spectrum of the evolved 
perturbations can be roughly separated into five regimes, having initially power-law shape 
and a cutoff scale as shown in Figure 1. 

According to these results, we can first separate the power spectrum roughly into three 
regimes on the scales larger than the cutoff scale: one is the linear regime (k < k n f. Regime 
1). The value of the power index in this regime remains the initial power index, n. The 
second regime is the single-caustic regime (k n i < k < k sn i: Regime 2). In this regime, 
the power index becomes —1 and is independent of the initial conditions. This result is 
caused by the appearance of caustics at this scale, and these caustics determine the power 
index of the power spectrum in this regime. The third regime is the multi-caustics regime 
(k sn i < k < k cut : Regime 3). In this regime, caustics on various scales of waves determine 
the power spectrum. Therefore, on these scales, the power spectral index ji can be traced 
back to initial conditions. 

On the scales smaller than the cutoff scale the problem is more involved. Nevertheless, 
we can roughly separate the power spectrum into two more regimes: after the appearance 
of the first caustic, the power index in the regime of the scale smaller than the cutoff scale 
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becomes —1. This value is, as mentioned above, caused by the appearance of caustics. How- 
ever, thereafter, more and more caustics appear and so the separation of caustics becomes 
smaller and smaller. Even after the appearance of many caustics, on scales smaller than the 
characteristic separation of the caustics (k > k cs : Regime 5), the power index of the power 
spectrum can be obtained by the density profile around the singular point: the power index 
on these scales also becomes —1, a result that can be derived from catastrophe theory (or the 
Lagrange-singularity theory, respectively), see, e.g., Arnol'd et al. 1982, Kotok & Shandarin 
(1988), Gouda & Nakamura (1988, 1989), Yano & Gouda (1998). On the other hand, on 
the scales larger than the characteristic separation of the caustics, in the fully developed 
'virialized' regime {k cut < k < k cs : Regime 4), the distribution of the singular points deter- 
mines the power index on these scales instead of the density profile around one singularity: 
we may say that the smoothed density profile with the smoothing scale (k cut < k < k cs ) 
determines the power index. For a comprehensive description of the above discussion, see 
Yano & Gouda (1998). 

In this paper we focus on the fully developed 'virialized' multi-stream regime (k cut < 
k < k cs ) to clarify the physical process of the appearances of the caustics, and we are going 
to estimate a value for the power index of the power spectrum in this regime. For this 
purpose we study the time-intervals of the appearances of the caustics as well as the spatial 
distribution of the caustics at subsequent times. The initial condition is chosen as 'single- 
wave' perturbation defined by v = vq sin x, where v is the velocity, x is the position of a sheet, 
and Vq is an arbitrary constant. The evolution of such a perturbation determines the power 
spectral index in the 'virialized' regime, essentially Region 4 in Figure 1. Accordingly we 
succeed to find some universal rules about the time when each caustic appears, the spatial 
distribution of the caustics at subsequent times, and the distribution of mass for each caustic. 

Our paper is organized as follows. In §2, we will briefly explain the numerical method 
to calculate the time evolution of density perturbations in the expanding universe model. In 
§3 we show the numerical results providing some rules for the time-appearances of caustics, 
the spatial distribution of caustics at subsequent times, and the mass distribution for each 
caustic. Here, we also evaluate the power index of the power spectrum in the 'virialized' 
regime. The final section is devoted to a summary and a brief outlook. 

2. MODEL AND NUMERICAL METHOD 

We investigate the time evolution of a plane-symmetric density perturbation in an 
Einstein-de Sitter universe model by using a semi-numerical method. In the one-dimensional 
system many plane-parallel sheets move only in perpendicular direction to the surface of 
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these sheets. When two sheets cross, they are allowed to pass through eachother freely. In 
this sheet system, there is an exact solution until two sheets cross over as follows (Zel'dovich 
1970, Sunyaev k Zel'dovich 1972, Doroshkevich et.al. 1973, Zentsova k Chernin 1980, 
Buchert 1989, Yano k Gouda 1998): 

x = q + B 1 (t)S 1 (q) + B 2 (t)S 2 (q), 

v = B 1 (t)S l (q) + B 2 (t)S 2 (q) , (1) 

where q and x are the Lagrangian and the (comoving) Eulerian coordinates, respectively. 
Here, Si(q), and S 2 (q) are arbitrary functions of q. Bi(t), and B 2 (t) are the growing mode 
and the decaying mode of linear perturbation solutions, respectively. We are considering 
the Einstein-de Sitter background universe model, i.e., Bi(t) = a and B 2 {t) = cTa, where 
a is the scale factor. The velocity v is the peculiar-velocity normalized by the scale-factor. 
We can compute the crossing time of all neighboring pairs of sheets. We use the shortest of 
these crossing times as a time step. Then, we can compute the new positions and velocities 
for all sheets at this crossing time. After two sheets cross, we exchange the velocities of two 
sheets that just crossed. Then we obtain again Sx(q), S 2 (q), and therefore exact solutions as 
follows: 

3 2 
Si(q) = -cT 1 ^ - q) + 

2 3 2 s 

S 2 (q) = -a2(x - q) - -a^a^v . (2) 

O 

These new exact solutions can be used until again two sheets cross. In this way we obtain 
the exact loci of the sheets by coupling these solutions. This semi-numerical method has 
good accuracy, because we connect exact solutions. 

In our calculation we are going to use 2 12 sheets. Periodic boundary conditions are fixed 
at a length of 2n. Then we consider the evolution of the 'single-wave' perturbation whose 
wave number is 1. 



3. RESULTS 

In this section we calculate the time-intervals for the appearances of caustics, the spatial 
distribution of caustics at subsequent times, and the mass distribution for each caustic. 
Accordingly we find some rules for the evolution of the caustic distribution. Using these 
rules, we are going to estimate the power index of the power spectrum in the 'virialized' 
regime. 
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3.1. Universality rules and constants for the occurence of caustics 

We examine the time evolution of the 'single-wave' perturbation setting Si(q) = sing 
and S2(q) = in Equations (1). Then, comoving positions x and scaled peculiar-velocities 
v of each sheet are simplified to: 

x = g + Bi(i)sin(g) , 

v = Si(*)sin(g) , (3) 

see Figure 2-a. Hereafter we call this case the single-wave case. 

The first caustic has just appeared at the center of the x-axis (x = 0) in Figure 2-b. 
After the first appearance of a caustic, more caustics appear at the center (phase mixing 
starts in phase space) (Figures 2-c - 2-d) (Doroshkevich et. al. 1980; Melott 1983; Gouda 
& Nakamura 1989; Yano & Gouda 1998). According to this cascade of appearing caustics, 
the density distribution evolves as in Figures (3-a - 3-d). Figures 2-a, 2-b, 2-c, and 2-d 
correspond to Figures 3-a, 3-b, 3-c, and 3-d, respectively. As we see from Figures 3-a, 3-b, 
3-c, and 3-d, a hierarchy of nested density peaks appears in the course of time. 

In order to discover rules for the distribution of the caustics at an arbitrary time, we 
sort the time-intervals in the appearances of the caustics, the spatial distribution of the 
caustics, and the mass distribution for each caustic at a certain time as follows. 

First, we examine the time-intervals related to the appearances of the caustics. The 
solid line in Figure 4 shows the time when caustics appear at the center of the x-axis (a so- 
called A3-type singularity). The dashed line shows the time-intervals for the appearances of 
caustics. As we can see from Figure 4, these time-intervals are almost constant. In order to 
clarify this point, we introduce the ratio of those intervals defined by Tj = (t i+2 — U+i)/ {ti+i — 
ti), where tj is the time when the i-th caustics appears. Then, the relation between Tj and 
the time when the i-th caustic appears is plotted in Figure 5. We see a loss of accuracy 
after around the 10th caustic appeared. Nevertheless, independent of the number of sheets, 
we see that the value Tj quickly converges to I. We can understand why the ratio becomes 
1 as follows: according to the time evolution of the density distribution, more and more 
caustics appear at the center of the x-axis and these caustics may be viewed as forming 
a bound cluster. After the evolution of this cluster decouples from the expansion of the 
universe model, sheets around the center of the cluster behave like a harmonic oscillator 
due to the fact that an enhanced gravitational potential forces particles with velocities less 
than the escape velocity from the cluster to rebound. The periods of the oscillation of the 
sheets finally becomes independent of the amplitude. Therefore, it is reasonable that the 
time-intervals of the appearances of caustics becomes constant, i.e., the value of Tj becomes 
1. 



- 6- 



Second, we consider the spatial distribution of caustics. In Figure 6 we plot the spatial 
positions of caustics when the N-th caustic appears at the center of the x-axis. Here we 
note that the number of appearances of caustics can be used instead of time, because we 
have already shown that the intervals may be considered practically constant. We define 
the position of a caustic, X(i,j), as the i-th position of a caustic at the time when the 
j-th caustic appears. Then, for example, X(i,i) represents the position of the i-th caustic 
that has just appeared as A3-type at the center of the x-axis. Now, we examine the spatial 
positions of caustics at a certain time (N = N ). We adopt N = 10 in our analysis because 
the accuracy of the positions of caustics decreases for N > 10 as shown in Figure 5. The 
solid and dashed lines in Figure 7 denote the spatial position of the N-th caustic, X(N, 10), 
and the spatial interval between the (N-l)-th and the N-th caustics, X(N—1, 10)—X(N, 10), 
respectively. Furthermore we depict the corresponding ratios: the solid and dashed lines in 
Figure 8 depict the ratio of spatial positions X(N, 10)/X(N — 1, 10) and the ratio of spatial 
intervals (X(N - 1, 10) - X(N, 10))/(X(N - 2, 10) - X(N - 1, 10)), respectively. In Figure 
8 we see no constant in the spatial positions or spatial intervals of caustics at first glance. 

However, we see a constant when we plot positions of the nearest caustic from the 
center, X(N - 1, N), or the ratio X(N - 1, N)/X(N - 2, N - 1). In Figure 9, the relation 
between the position X(N — 1, N) and the time when the N-th caustic appears is plotted. 
The ratio of positions X(N - 1, N)/X(N - 2, TV - 1) is plotted in Figure 10. In Figure 10, 
Ri represents the ratio between the position of the nearest caustic from the center and that 
of the second nearest caustic from the center X(N — 1, N)/X(N — 2, N). R 2 represents the 
ratio between the nearest caustic from the center at the (N-l)th appearance of a caustic and 
that at the Nth appearance of a caustic X(N — 1, N)/X(N — 2, N — 1). Both ratios feature 
constants: the ratio R\ and R 2 have values of around 0.15 and 0.3, respectively. 

If we plot positions of caustics by using a logarithmic scale of spatial positions and a 
linear scale of time, we can see some "lattice" as shown in Figure 11. This apparent regularity 
is the result of the existence of constants in the ratios. Here we note that calculated positions 
are not exact for smaller than about 10~ 6 because of finite sheet resolution. Therefore, the 
"lattice" is destroyed after around the 10-th appearance of caustics. 

Third, we investigate the mass distribution for each caustic at a certain time (N = N ). 
Here, we define the mass, M(i, N), as proportional to the number of sheets which are located 
inbetween X(i,N) and X(i + 1,N) at the time, when the Nth caustic appears. The schematic 
picture of M is shown in Figure 12. 

Then, we obtain the mass distribution for each caustic at N = 10, which is shown in 
Table 1 and Figure 13. We deduce the following relation between % and M for large values 



-7- 



of v. 

M(i, N) \ I!' XI . where R M ~ 0.5 . (4) 

Since we know the spatial distribution of the caustics and the mass distribution for the 
caustics, we can construct the smoothed density profile to calculate the power index of the 
power spectrum in the 'virialized' regime. 



3.2. Evaluation of the power spectrum in the 'virialized' regime 

In this subsection we estimate the index of the power spectrum in the 'virialized' regime, 
using the rules which we have found in the previous subsection. Assuming that the spatial 
distribution of caustics obeys the law of a geometric series, 

X(i, N ) = a a(N y , (5) 

that is, 

lnX(i, N ) - lna , r , 

* = i , 6 

ma 

and also that the mass distribution for each caustic obeys the law of another geometric series, 

M(i,N ) = b b(N ) i , (7) 

we find the density profile as: 

n(X) = M(hNo) - b ° (b_Y lna ° /lna x l^M (R) 

Q[ ' ~ X(i,N Q )-X(i + l,N Q ) (l-a)a \a) ' [) 

Its Fourier transform gives 

5 k oc J Q(X)e lkX dX 
oc Jx^e ikX dX 

oc AT 1 - 1 ^ . (9) 
Finally we obtain the power spectrum as: 

P(k) = |4| 2 oc k- 2 - 2 ^ . (10) 



In fact, we have found a = R\ = 0.15 and b = Rm = 0.5 from the previous subsection, so 
that we obtain the power spectrum in the explicit form: 

P(k) oc k\ where v = -0.73. (11) 
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On scales smaller than the characteristic scale k > k cs (Regime 5 in Figure 1), the power 
index of the power spectrum becomes —1. This is because the density profile around one 
singular point determines the power spectrum on this scale. 

On scales larger than the characteristic scale k < k cs , on the other hand, the density 
profile in the spatial interval, where the singular points lie closer than the scale of this 
interval, is smoothed, and the power index of the power spectrum deviates from —1. The 
larger scales we consider, the more the density profile is smoothed. The power index in the 
'virialized' regime k cut < k < k cs , accordingly, is determined by the smoothed density profile 
(8), and therefore the result for the power index of the power spectrum under the smoothed 
density profile is given by (11). 

Our result shows that the power index in the 'virialized' regime k cut < k < k cs may be 
approximated by the index value v = —0.73. In addition, let us note that the wavenumbers 
k cs , which represent the characteristic separations of caustics, becomes larger and larger as 
time increases, since the characteristic separation of caustics becomes shorter and shorter. 

4. SUMMARY AND OUTLOOK 

We have investigated the evolution of density perturbations after the first appearance 
of a caustic, exemplified for a plane-symmetric 'single-wave' perturbation in an expanding 
universe model. Following Yano & Gouda (1998), who showed that the power spectrum can 
be roughly separated into five regimes with different power-law index, we have explored the 
fourth regime. The detailed analysis of a many-sheet system is needed in this regime in 
order to determine the time-evolution and the asymptotics of the power spectral index. 

We have argued that the system of caustics displays some universality, after some time, 
with respect to regularity of time-intervals in the appearances of caustics, and the spatial 
distribution of caustics at subsequent times. 

First, we have noticed the existence of a rule for the time-intervals of the appearances 
of caustics: the ratio of intervals quickly approaches 1, a fact that we have made plausible. 
Second, we have found some rules hidden in the spatial distribution of caustics: the ratio 
between a position of the nearest caustic from the center and that of the second nearest 
caustic from the center Ri = X(N — 1, N)/X(N — 2, N) approaches a value of around 0.15. 
Furthermore, the ratio between the nearest caustic from the center at the (N-l)th appearance 
of a caustic and that at the Nth appearance of a caustic R 2 = X(N — 1,N)/X(N — 2,N — 1) 
attains a value of around 0.3. 
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Third, we have found a rule for the mass distribution for each caustic, that is, we have 
derived the density distribution of a cluster. The mass distribution for each caustic can be 
fit to a line with the geometric progression of ratio 0.5. 

Furthermore, using the above constants, we have found the smoothed density profile, 
and finally have evaluated the power index of the power spectrum in the 'virialized' regime 
to be —0.73. This value is different from —1 that is determined by the density profile around 
one singular point. In other words, the power spectum in the region below the cutoff scale 
can be divided into two clearly distinct regions, the 'single-caustic' regime and the 'virialized' 
regime. These results are consistent with the numerical results by Yano & Gouda (1998). 

The discovery presented in this research note, namely that there exist remarkable rules 
for the distribution of caustics, points towards several interesting theoretical implications for 
long-term properties of self-gravitating systems. As in studies of dynamical systems we here 
deal with a hierarchy of bifurcations. A fruitful way to explore universality properties further 
is to use renormalization group techniques similar to those developed for the Feigenbaum 
scenario of pitchfork bifurcations (see: Feigenbaum 1978, 1979; Derrida et al. 1979; Collet 
& Eckmann 1980). With this note we have just provided a hint to the existence of universal 
constants. A detailed renormalization group and fixpoint analysis is needed to support these 
findings. In this context an interesting question has to remain open: the existence of chaotic 
regimes. 

Besides this theoretical interest we think that our results may be useful to check the 
accuracy of numerical simulations concerning their stability under long-time integration. 

In this paper we have focused on one non-trivial regime, the 'virialized' regime (Regime 
4 in Figure 1), and have succeeded to evaluate the power index of the power spectrum in 
this regime semi-numerically. In order to complete the understanding of the evolution of the 
power spectra of density fluctuations, we have also to understand the behavior of the power 
index in the multi-wavenumbers regime (Regime 3 in Figure 1). 

Furthermore, we could employ the systematic framework of the velocity moment hierar- 
chy to eventually derive analytically universal behavior of the density profile; also, nonlinear 
statistical methods beyond the power spectrum would give more refined information. 
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Fig. 1. — Schematic representation of the evolved power spectrum for a general initial power 
spectrum after the first appearance of caustics. 



Fig. 2. — (a) Distribution in phase space for the 'single-wave' case at the initial time (a), 
at the first appearance of caustics (b), at the second appearance of caustics (c), and at the 
fifth appearance of caustics (d). 



Table 1: Mass distribution for each caustic when the 10th caustics appear. 



order of caustics 


position of caustics 


numbers of sheets 


1 


0.2736580 


1094 


2 


0.1288644 


2021 


3 


0.0507813 


1888 


4 


0.0170624 


1289 


5 


0.0051794 


763 


6 


0.0014272 


414 


7 


0.0003467 


204 


8 


0.0000683 


95 


9 


0.0000081 


42 
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Fig. 3. — (a)Density distribution for the 'single-wave' case at initial time (a), and at the 
first (b), the second (c) and the fifth (d) appearance of caustics. 




Fig. 4. — Times and intervals for the appearance of caustics at the center of the x-axis. The 
solid and dotted lines represent times and intervals, respectively. 
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Fig. 5. — Ratio of intervals in the appearance of caustics. The solid, long-dashed, short- 
dashed, dotted, and dash-dotted lines are cases for different numbers of sheets: n=4095, 
n=2047, n=1023, n=511, and n=255, respectively. 
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Fig. 6. — Spatial distribution of caustics at the N-th appearance of caustics (the bifurcation 
diagram) . 
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Fig. 7. — Spatial positions and intervals of caustics at the 10-th appearance of caustics. The 
solid and dashed lines represent the spatial positions and intervals, respectively. Smaller N 
represents a position that is further away from the center. 
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Fig. 8. — Ratio of spatial positions and intervals of caustics at the 10-th appearance of caus- 
tics. The solid and dashed lines denote ratios of spatial positions and intervals, respectively. 
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9. — Spatial position of the caustic that is nearest to the center of the x-axis X(N— 1, N). 
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Fig. 10. — Ratio of positions of caustics. R\ represents the ratio between the position of 
the nearest caustic from the center and that of the second nearest caustic from the center 
X(N — 1,N)/X(N — 2,7V). R 2 represents the ratio between the nearest caustic from the 
center at the (N-l)th and that at the Nth caustic appearance, X(N-1, N)/X(N-2, N- 1). 
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Fig. 11. — Spatial positions of caustics at the N-th appearance of caustics. 
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Fig. 12. — Schematic picture of the density distribution and the definition of the mass for 
the i-th caustic (the shaded region). 
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Fig. 13. — Mass distribution contributing to each caustic at the 10-th appearance of caustics. 



